function [S,F,obj,obj2] = MyFunction(Z_ini,A,num_cluster,lambda1,lambda2,theta,rho,max_iter,delta,alpha,beta)

%--------------------------Initialize---------------------------%
%---------initialize S,B------------------
S = Z_ini;
B = S;
%---------initialize J--------------------
for iv = 1:length(S)
    J{iv} = zeros(size(S{iv}));
end
%-------------------------------------------
Nsamp = size(S{1},1);
nv = length(S);

S_tensor = cat(3, S{:,:});
d = size(S_tensor);
dtemp=d(1)*d(2)*d(3);
j=zeros(dtemp,1);
b=zeros(dtemp,1);

miu = ones(nv,1) / nv;

directions = [1];
n = length(directions);

sX = zeros(Nsamp);
knn_idx = true(Nsamp);
up_knn_idx = triu(knn_idx);
jiaquan = ones(nv,1);
zr = 10e-10;
rho2 = 0.0001;

%-----------initialize G,Gamma-----------------%
for i = 1:n
    index        = directions(i);
    G{index}     = porder_diff(S_tensor,index); 
    Gamma{index} = zeros(d); 
end

%---------------initialize U-----------------%
U = zeros(Nsamp);
for v = 1:nv
    U = U + S{v};  
end
U = U/nv;
for j = 1:Nsamp
    if sum(U(j,:)) ~= 0
        U(j,:) = U(j,:)/sum(U(j,:));
    end
    
end;
sU = (U+U')/2;
D = diag(sum(sU));
L = D - sU;
[F, ~, evs] = eig1(L, num_cluster, 0);
U = sU;



for iter = 1:max_iter
    S_pre = S;
    B_pre = B;
    F_pre = F;

        % ------------------- C(i) -------------------- %
        W = alpha*ones(nv) - diag(alpha*ones(1,nv)) + diag(beta*ones(1,nv));
        M = W.*(miu*miu') + diag(miu);

        commom_baA = zeros(Nsamp);
        for v = 1:nv 
            baA{v} = alpha*miu(v)*S{v};
            special_baA{v} = beta*miu(v)*S{v};
            commom_baA = commom_baA + baA{v};
        end
        for v = 1:nv 
            true_baA{v} = commom_baA - baA{v} + special_baA{v};
            temp = full(miu(v)*(U + true_baA{v})); % 
            P{v} = temp(up_knn_idx);
        end
        
        right_p = cat(2, P{:})';
        if det(M) == 0
            solution = (pinv(M) * right_p)';
            fprintf('------------')
        else
            solution = (M \ right_p)';
        end
        solution(solution<0) = 0;
        for v = 1:nv
            temp = solution(:,v);
            C{v} = zeros(Nsamp);
            C{v}(up_knn_idx) = temp;
            C{v} = max(C{v}, C{v}');
            C{v} = min(C{v}, S{v});
            
        end
        % ------------------- S(i) -------------------- %
        commom_baE = zeros(Nsamp);
        for i = 1:nv  
            baE{i} = alpha*miu(i)*(S{i}-C{i});
            special_baE{i} = beta*miu(i)*(S{i}-C{i});
            commom_baE = commom_baE + baE{i};
        end
        for i = 1:nv 
            true_baE{i} = commom_baE - baE{i} + special_baE{i};
            tem = full(miu(i)*true_baE{i});
            zeta{i} = tem;
            zeta{i} = max(tem, tem');
        end

        for iv = 1:nv
            linshi_S = (Z_ini{iv}.*A{iv}+theta*0.5*B{iv}-0.5*J{iv}-0.5*zeta{iv})./(0.5*theta+A{iv});
            Z1 = zeros(size(linshi_S));
            for is = 1:Nsamp
                ind_c = 1:Nsamp;
                ind_c(is) = [];
                Z1(is,ind_c) = EProjSimplex_new(linshi_S(is,ind_c));
            end
            S{iv} = Z1;
        end
        clear Z1 linshi_S linshi_M
        % ---------------------------- B -------------------------- %
        S_tensor = cat(3, S{:,:});
        %% FFT setting
        T = zeros(d);
        for i = 1:n
            Eny = diff_element(d,directions(i));
            T   = T + Eny; 
        end
        %% Update B -- solve TV by FFT 
        HB = zeros(d);
        for i = 1:n
            index = directions(i);
            HB = HB + porder_diff_T(theta*G{index}-Gamma{index},index); 
        end
        J_tensor = cat(3, J{:,:});
        B_tensor = real( ifftn( fftn( theta*S_tensor+J_tensor+HB)./(theta*(1+T)) ) );

        % ------------------------- Gi --------------------------------- %
        for i = 1:n
            index = directions(i);
            linshi_add=porder_diff(B_tensor,index)+Gamma{index}/theta;
            [G{index},tnn_G{index}] = solve_G(linshi_add,n*theta/lambda1,delta);
        end
       % ------------------------- J{iv} --------------------------------- %
        s=S_tensor(:);
        b=B_tensor(:);
        j=j+theta*(s-b);
        for iv = 1:nv
            B{iv} = B_tensor(:,:,iv);
            J_tensor=reshape(j,d);
            J{iv} = J_tensor(:,:,iv);
        end

        % --------------------------- F ---------------------------------- %
        U1 = U;
        U1 = (U1 + U1')/2;
        D = diag(sum(U1));
        L = D - U1;
        F_old = F; 
        [F, ~, ev]=eig1(L, num_cluster, 0);

       % ---------------------------- U ---------------------------------- %
        U = updateU(C, F, miu, lambda2, Nsamp, nv);

       % ------------ Update mulipliers:  Gamma, theta, miu----------- %
        for i = 1:n
            index = directions(i); 
            Gamma{index} = Gamma{index}+theta*(porder_diff(B_tensor,index)-G{index});
        end

        theta = min(theta*rho, 1e10);
    
        for v = 1:nv
            distUC = norm(U - C{v},'fro')^2;
            if distUC == 0
                distUC = eps;
            end;
            miu(v) = 0.5/sqrt(distUC);
        end
        clear B_tensor S_tensor J_tensor

       % ------------------------------ obj ------------------------------ %

        diff_S = 0;
        diff_B = 0;
        linshi_obj  = 0;
        %% check convergence
        for iv = 1:nv
            leq{iv} = S{iv}-B{iv};
            diff_S = max(diff_S,max(abs(S{iv}(:)-S_pre{iv}(:))));
            diff_B = max(diff_B,max(abs(B{iv}(:)-B_pre{iv}(:)))); 
        end
        leqm = cat(3, leq{:,:});
        leqm2 = max(abs(leqm(:)));
        clear leqm leq
    
        err = max([leqm2,diff_S,diff_B]);
        % fprintf('iter = %d, theta = %.3f, err = %.8f, difS = %.8f, diffSB = %.8f\n'...
        %     , iter,theta,err,diff_S,leqm2);
        obj(iter) = err; 
        obj2(iter) = leqm2;    
        if iter>5 & leqm2 < 1e-6
%             iter
            break;
        end    
    
    end

function U = updateU(C, F, mu, lambda, numSmp, numView)
    dist = L2_distance_1(F',F');
    U = zeros(numSmp);
    for i=1:numSmp
        c0 = zeros(1,numSmp);
        for v = 1:numView
            temp = C{v};
            c0 = c0 + mu(v)*temp(i,:);
        end     
        idxa0 = find(c0>0);
        ci = c0(idxa0);
        ui = dist(i,idxa0);
        cu = (ci - 0.5*lambda*ui)/sum(mu);
        U(i,idxa0) = EProjSimplex_new(cu);
    end;
end



























end

